Osculating orbits in Schwarzschild spacetime, with an application to extreme 

mass-ratio inspirals 
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We present a method to integrate tiie equations of motion tiiat govern bound, accelerated orbits 
in Schwarzschild spacetime. At each instant the true worldline is assumed to lie tangent to a 
reference geodesic, called an osculating orbit, such that the worldline evolves smoothly from one 
such geodesic to the next. Because a geodesic is uniquely identified by a set of constant orbital 
elements, the transition between osculating orbits corresponds to an evolution of the elements. In 
this paper we derive the evolution equations for a convenient set of orbital elements, assuming 
that the force acts only within the orbital plane; this is the only restriction that we impose on the 
formalism, and we do not assume that the force must be small. As an application of our method, we 
analyze the relative motion of two massive bodies, assuming that one body is much smaller than the 
other. Using the hybrid Schwarzschild/post-Newtonian equations of motion formulated by Kidder, 
Will, and Wiseman, we treat the unperturbed motion as geodesic in a Schwarzschild spacetime with 
a mass parameter equal to the system's total mass. The force then consists of terms that depend on 
the system's reduced mass. We highlight the importance of conservative terms in this force, which 
cause significant long-term changes in the time-dependence and phase of the relative orbit. From 
our results we infer some general limitations of the radiative approximation to the gravitational 
self-force, which uses only the dissipative terms in the force. 

PACS numbers: 04.20.-q, 04.25.-g, 04.25. Nx, 04.40.-b 



I. INTRODUCTION 

A. Orbital motion in curved spacetime 

Analysis of accelerated orbits in curved spacetime has 
historically focused on the post-Newtonian regime (see 
Ref. [ESQ for general reviews of the post-Newtonian 
formalism), since observations of orbital motion have his- 
torically been limited to weak-field systems such as the 
solar neighborhood and binary pulsars. However, the 
advent of gravitational-wave astronomy has recently ne- 
cessitated an analysis of accelerated orbits in strongly- 
curved spacetimes. The primary examples of such or- 
bits are extreme mass-ratio inspirals (EMRIs), in which 
a small compact body of mass m spirals into a supermas- 
sive black hole of mass M ^ m. Such systems promise to 
be excellent sources of gravitational waves for the space- 
based detector LISA However, accurate predictions 
of the emitted waveforms must account for the effect of 
the compact body's gravitational field on its own mo- 
tion. The compact body induces a metric perturbation 
hai3 = {m/M)h^^^p + 0(m/M)2. Although the motion 
of the particle may be described as a geodesic in the 
perturbed spacetime, it is more simply treated as an ac- 
celerated worldline in the background spacetime of the 
unperturbed black hole. The cause of the acceleration 
is thus interpreted as a gravitational self-force derived 
from a regularized form of the field hafj. This force was 
first formally calculated to first order in m/M by Mino, 
Sasaki, and Tanaka Q, and later by Quinn and Wald 6] 
(see Ref. Q for a review of recent developments). Other 
possible effects on the inspiraling particle, such as tidal 
perturbations of the central black hole, spin-orbit and 



spin-spin couplings, electromagnetic interactions, and so 
on, can also be treated as forces acting on the body. 

Although significant progress has been made in cal- 
culating these effects (see Ref. Q for a recent review of 
work on EMRIs), there has been no attempt to formulate 
a general method of determining and characterizing the 
resulting motion. Implementing the first-order gravita- 
tional self-force brings a particular difficulty: The self- 
force on a particle is a functional of the particle's world- 
line, which for the first-order calculation is assumed to be 
a geodesic. However, the true motion is never geodesic, 
because of the self-force. Thus, the effect of the self-force 
must somehow be determined with reference to a ficti- 
tious geodesic worldline. 

In this paper we present a method to integrate the 
equations of motion that govern accelerated motion in 
Schwarzschild spacetime. The method can be used for a 
wide class of perturbing forces; the only restrictions are 
that the force must keep the orbital motion bounded be- 
tween a minimum and a maximum radius (the method is 
not suitable for the final portion of an orbit that plunges 
into the black hole), and that it must be acting within the 
plane of the orbit (although the method could be easily 
extended to accommodate non-planar motion). Within 
these restrictions the force is arbitrary, and in particular, 
it is not assumed to be small. Our method is a rela- 
tivistic extension of the traditional method of osculating 
orbits, also called the method of variation of constants, 
in Newtonian celestial mechanics (see, e.g., Refs. [ol. [loj). 
In this method the true worldline z{X) is taken to lie tan- 
gent to a geodesic zg{X) at each value of the orbital pa- 
rameter A, such that the true orbit moves smoothly from 
one geodesic to the next. The instantaneously tangential 
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geodesies are referred to as osculating orbits (meaning 
"kissing orbits"). A geodesic is characterized by a set 
of constants I^, called orbital elements, and the tran- 
sition between osculating orbits corresponds to changes 
in these elements; thus, the method of osculating orbits 
amounts to parametrizing the true worldlinc as an evolv- 
ing geodesic with dynamical orbital elements /"^(A). 

Because it explicitly determines the position and veloc- 
ity of a tangential geodesic at each instant, this method 
explicitly provides the information necessary to calcu- 
late the first-order gravitational self- force at each instant. 
Our method is therefore very well suited to the gravita- 
tional self- force problem. It also has several more general 
advantages. First, because the orbital elements are con- 
stant on a geodesic, the method clearly separates pertur- 
bative from non-perturbative effects. (Throughout this 
paper the accelerated motion of the particle is referred to 
as a perturbation of the geodesic motion. However, this is 
only to distinguish effects of acceleration from effects on 
a geodesic; the "perturbation" need not be small.) Sec- 
ond, although the orbital elements are equivalent to the 
set of initial conditions, they are typically chosen so as 
to provide direct geometric information about the orbit. 
If the perturbing force is very weak, then the perturbed 
orbit will he very close to a geodesic for a long period of 
time, and changes in the orbital elements will character- 
ize changes in the geometry of the orbit. Thus, although 
our method is exact, it is perhaps most useful in the con- 
text of small perturbations. Third, the orbital elements 
divide into two classes. The first class includes the prin- 
cipal orbital elements; these are equivalent to constants 
of the motion such as energy and angular momentum, 
and they determine the geodesic on which the particle is 
moving. The second class includes the positional orbital 
elements, which determine the particle's initial position 
on the selected geodesic, as well as the geodesic's spa- 
tial orientation. Generally speaking, long-term changes 
in the principal orbital elements are produced by dis- 
sipative terms in the perturbing force, while long-term 
changes in the positional elements are produced by con- 
servative terms. Thus, this division into two classes al- 
lows one to easily separate conservative from dissipative 
effects of the perturbing force. 

We note that this general idea of characterizing or- 
bital evolutions by changes in the "constants" of motion 
has been used frequently in analyzing the effects of ra- 
diation reaction. Such analyses have typically focused 
on changes in the principal elements alone, neglecting 
the changes in positional elements, and rarely mention- 
ing the general framework of osculating orbits. However, 
there have been at least two notable generalizations of the 
method of osculating orbits from Newtonian to relativis- 
tic mechanics: the adaptation of the method by Damour 
et al. to post-Newtonian binary systems [l^], and 
the formulation proposed by Mino for orbits around a 
Kerr black hole The formulation by Damour et al. 
is complete and easy to implement, but it is limited to 
the post-Newtonian regime. Mino's formulation is valid 



for arbitrary bound orbits in Kerr, and it was undoubt- 
edly useful for Mino's own purposes, but we believe that 
a concrete implementation of his method would not be 
very practical. The reason is that Mino expresses the 
orbits as formal Fourier expansions with unknown con- 
vergence behavior, in terms of coefficients that would be 
difficult to calculate in practice. It may well be that 
the complexity of geodesies in Kerr make a more practi- 
cal parametrization impossible, but as we shall demon- 
strate in this paper, we can do much better for orbits in 
Schwarzschild spacetime. Given the limitations of pre- 
vious work, we believe that it is timely to present here 
a practical formulation of the method of osculating or- 
bits for bound motion in Schwarzschild spacetime. We 
shall first present an outline of the general method in rel- 
ativistic mechanics and its connection to the traditional 
method in Newtonian mechanics, and we shall next spe- 
cialize the method to the case of bound orbital motion in 
Schwarzschild spacetime. 

B. Test case: post-Newtonian binaries 

We demonstrate the usefulness of our method by ap- 
plying it to the relatively simple system of two compact 
bodies of mass mi and m2 3> rui in the post-Newtonian 
regime. The equations of motion for the spatial positions 
Xi and of the bodies have been determined in har- 
monic coordinates to 3.5PN order (i.e. of order {v/cf 
beyond the Newtonian description) [31. Conservative 
terms appear at IPN, 2PN, and 3PN orders, and dis- 
sipative terms appear at 2.5PN and 3.5PN orders. Since 
the essential features of the problem are already present 
at 2.5PN order, we truncate the equations at that order 
for simplicity. These equations are valid for arbitrary 
mass ratios, but we focus on the extreme case in order to 
link our results to the self-force problem. 

In order to analyze this system of equations with our 
method of osculating orbits, we use the hybrid equa- 
tions of motion constructed by Kidder, Will, and Wise- 
man These equations take the schematic form 

^ = -^(l + SCHW + /iPF). (1) 

The spatial separation vector S connects the 

two bodies, and M = mi -\- m2 and /i = mim2/M are re- 
spectively the total mass and reduced mass of the system. 
The terms in SCHW are the exact relativistic corrections 
to Newton's law in a Schwarzschild spacetime of mass M, 
so that = — ^ (1 -I- SCHw) is the exact equation for 
a test particle in that spacetime. The terms in /iPF are 
those in the post-Newtonian expansion that depend ex- 
plicitly on the reduced mass of the system (pf stands for 
"perturbing force"). Since the extra terms introduced 
within SCHW are of 3PN order and higher, the hybrid 
equations remain correct at 2.5PN order. However, they 
differ from the usual post-Newtonian equations in that 
they become exact in the test- mass limit fj, 0. This 
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allows us to apply our method to the post-Newtonian 
system by taking our osculating orbits to be geodesies in 
the fictitious Schwarzschild spacetime of mass M, and by 
deriving our perturbing force from ^pf. 

The force derived in this way is a form of the grav- 
itational self-force, since it is produced by finite-mass 
effects. However, it differs nontrivially from the post- 
Newtonian limit of the relativistic self-force: First, the 
self-force is a gauge- dependent quantity which is typically 
calculated in the Lorenz gauge, while the hybrid equa- 
tions of motion are derived within the harmonic gauge 
of post- Newtonian theory. Second, the Lorenz gauge en- 
sures that the coordinates of the small body are defined 
in relation to the system's center of mass [15|, while here 
we use coordinates relative to the large mass. And third, 
our geodesies are in a fictitious Schwarzschild spacetime 
of mass M — mi+ TO2 and not in the background space- 
time of the second body (of mass 7712). The last two 
differences could be easily removed by formulating an al- 
ternative set of hybrid equations, but the gauge difference 
cannot be easily dealt with. 

Given these differences, our method of osculating ele- 
ments is used in this paper primarily as a practical means 
to integrate the hybrid equations of motion. Neverthe- 
less, the perturbing force that we derive and the gravita- 
tional self-force share many essential features. In partic- 
ular, the self-force can be expected to have conservative 
terms at OPN (the Newtonian level), IPN, and 2PN or- 
ders, etc., and dissipative terms at 2.5PN (corresponding 
to quadrupole radiation) and 3.5PN orders, etc.; our per- 
turbing force has exactly the same features, except for the 
Newtonian correction, which is implicitly accounted for 
by working in terms of total and reduced masses. Thus, 
we can hope to draw some reasonable conclusions about 
the action of the gravitational self-force from our simpli- 
fied analysis. 

Our focus will be on detailing the limitations and ambi- 
guities of two approximation schemes, following our anal- 
ysis of the post-Newtonian electromagnetic self-force in 
Refs. [1^ The first scheme of interest lies within 

the broad class of adiabatic approximations, which rest 
on the assumption that the accelerated orbit deviates 
only "slowly" from the geodesic orbit. In particular, they 
commonly assume that any period of the motion is much 
shorter than the radiation-reaction timescale of the in- 
spiral, allowing one to eliminate irrelevant short-term os- 
cillations and keep only secular effects. Based on this 
assumption, an explicit implementation of such an ap- 
proximation will typically involve some type of averaging, 
either in the form of direct averaging of the equations of 
motion or via a two-timescale expansion. For clarity, we 
will refer to this averaging method, which is just a specific 
type of adiabatic approximation, as a secular approxima- 
tion. Using the hybrid equations of motion, we show 
in Sec. Ill B that the secular approximation introduces 
ambiguities in the choice of (a) initial conditions and (b) 
the variable to be averaged over. Our results suggest that 
different choices can significantly affect long-term behav- 



ior, and our conclusion is that while the idea of a secular 
approximation is attractive, the precise construction of 
one presents significant difficulties. 

We shall also examine the (pseudo-adiabatic) radiative 
approximation, which uses the radiative (half-retarded 
minus half-advanced) solution to the linearized Einstein 
equation. As shown by Mino [l^, the self- force calcu- 
lated from the radiative field approximately reproduces 
the long-term dissipative effects of the true self-force. 
Largely based on this result, it was believed that the 
radiative approximation would produce a valid adiabatic 
approximation to the true evolution. This notion has 
led to a confusing nomenclature in the literature, in 
which adiabatic and radiative approximations are treated 
synonymously. Since the radiative approximation intro- 
duces errors beyond those of an adiabatic approximation 
we find it misleading to identify the two. We 
insist here that the radiative approximation is logically 
distinct from the class of adiabatic approximations intro- 
duced in the preceding paragraph. 

Due to its simplicity, the radiative approximation 
has been utilized by several groups in analyzing EM- 
RIs [TqI , [20} . Unfortunately, the radiative self- force ne- 
glects all conservative effects of the true self-force. In 
the framework of osculating orbits, this translates into 
neglecting long-term changes in the positional orbital el- 
ements. (Although Mino has given prescriptions for find- 
ing these long-term changes using only the radiative self- 
force fisl. [Tpl . his prescriptions are highly ambiguous in 
practice [2l|.) As pointed out in Ref. ^0], the radiative 
approximation may have some utility despite this error, 
and in particular, it may be sufficient to generate tem- 
plates for the detection of a gravitational-wave signal. 
But it is unlikely that it will be sufficiently accurate for 
reliable parameter estimation. Because of the potential 
usefulness of the approximation, determining its limita- 
tions is quite important. In this paper we find that ne- 
glecting conservative effects leads to long-term errors in 
the phase and time-dependence of the orbit; this agrees 
with and extends our earlier results Il7| . The errors 
in the time dependence are of particular importance, as 
they apply even to the evolution of the principal orbital 
elements. 



C. Organization of this paper 

In Sec. Ill A] we introduce the general method of oscu- 
lating orbits. We then restrict our analysis in Sees. Ill Bl 
and III CI to bound planar orbits in Schwarzschild space- 
time. Section III Bl presents a parametrization of 
bound geodesies in terms of five orbital elements, and 
Sec. Ill CI uses the osculation condition to find evolu- 
tion equations for these orbital elements. In the sec- 
ond part of our paper, we apply our method to the hy- 
brid Schwarzschild/post-Newtonian equations of motion, 
which are presented in Sec. IIII Al The results of using 
a secular or radiative approximation are then displayed 
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and discussed in Sec. IIIIBI 



II. THE METHOD OF OSCULATING ORBITS 

A. The osculation condition 

We first consider the completely general situation of 
a point particle moving on an arbitrary worldline 2" (A) 
parametrized by A. We define the acceleration or 
force per unit mass, acting on the particle via the equa- 
tion of motion 



From Eq. ^ we have that = which implies 



r, 



(2) 



where an overdot indicates a derivative with respect to 
the proper time t on the worldline. The normalization 
condition z°'za = —1 implies the orthogonality condition 
f^Za — 0, which will be essential for later calculations. 
The relation between /" and the Newtonian perturbing 
force is discussed in Appendix [X\ 

Using the relations i" = ^A and = ^^'^ + ^^^ 
the equation of motion becomes 

dV^ dz^dz^^ fdr 

dX^ ^"^ dX dX ■' \ dX 



.ix/jl. (3) 



where k = —X/)?. The first term on the right-hand side 
is due to the force acting on the particle, while the second 
term is present whenever A is a non-affine parameter. 

Our goal is to transform the equation of motion ([3]) 
into evolution equations for a set of orbital elements I^. 
That is, we seek a transformation {z", i"} — > I^. Letting 
zg.(/'^,A) be a geodesic with orbital elements /"^, the 



osculation condition states the following: 



z"(A) 
dz"" 



^S(/<(A).A). 

^(/-(A),A), 



(4) 
(5) 



where the partial derivative in the second equation holds 
fixed. These two equations assert that at each value 
of A we can find a set of orbital elements I^{X) such that 
the geodesic with those elements has the same position 
and velocity as the accelerated orbit. We can freely make 
this assertion because the number of orbital elements is 
equal to the number of degrees of freedom on the orbit. 

As a consequence of the osculation condition, all re- 
lations that are obtained using only algebraic manipula- 
tions of coordinates and velocities on a geodesic are also 
valid on the true orbit. However, it is important to note 
that K is altered by the acceleration of the worldline, be- 
cause it involves second derivatives. Hence, an expression 
for k(A) that is valid on an osculating geodesic will not be 
valid on the tangential accelerated orbit. Nevertheless, 
A = for an affine parameter A on both orbits, so affine 
parameters remain affine. 

Now, combining the osculation condition with the 
equations of motion generates evolution equations for I^. 



dz" „ dz 

Comparing this result with Eq, 



a^g di^ 



where the index A is summed over. 
, we find 



dz^ dl^ 
dI^~dX 



0. 



Furthermore, Zq satisfies the geodesic equation 



9A2 



pa 



Ozq Ozq 



dX dX 



.c(A)^, 



(6) 



(7) 



where kg (A) is the measure of non-affinity of A on the 
geodesic. Subtracting this geodesic equation from the 
equation of motion ^ and using Eq. ^ to remove the 
Christoffel terms, we obtain 



dX^ 



9A2 



[«(A) 



-g(A)]^. (8) 



But differentiating Eq 

f d dz^\ di^ 
{W^ dX i d\ 



m yields ^ = ^ 

Comparing these results, we find 



d 



dz^ 



dl^ dX 



dl^ 
~dX 



KA)-«g(A)]^. 



(9) 



Equations © and Q form a closed system of first- 
order differential equations for the orbital elements I^. 
Two sources of change in the orbital elements are appar- 
ent: a direct source due to the perturbing force and 
an indirect source due to the change in the affinity of the 
parametrization of the accelerated orbit. Determining 
this second effect in practice may be somewhat difficult. 
However, if we use the affine parameter A = r then the 
equations simplify to 



dz"' 



G jA 



0, 

r. 



(10) 

(11) 



These equations can be easily inverted to solve for the 
derivatives /"^, which is done in Sec. IH CI If a non-affine 
parameter A is required in a specific application, one may 
easily find by multiplying the above equations by ^ , 
which will also be done in Sec. IH CI 



B. Geodesies in Schwarzschild spacetime 

We now focus on the specific case of bound or- 
bits in Schwarzschild spacetime. The osculating or- 
bits in this case are bound geodesies, for which we 
use the parametrization presented in the text by Chan- 
drasekhar and described in detail in Ref. "251. This 
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parametrization is given in Schwarzschild coordinates 
and can be easily derived as follows. 

Because of the spherical symmetry of the 
Schwarzschild spacetime, we can freely set 6 = 7r/2. 
The geodesic equations in a Schwarzschild spacetime 
with mass parameter M can be easily solved for the 
remaining coordinates to find 



E/F, 
L 



(12) 
(13) 

(14) 



where F — 1 — 2M/r, E and L are constants equal to 
energy and angular momentum per unit mass, respec- 
tively, the effective potential is UcS = F{1 + Ljr^), and 
an overdot represents a derivative with respect to the 
proper time r on the orbit. 

We are interested in bound orbits that oscillate be- 
tween a minimal radius r\ and a maximal radius ^2, re- 
spectively referred to as periapsis and apoapsis. Adapt- 
ing the tradition of celestial mechanics, we define the 
(dimensionless) semi-latus rectum p and the eccentricity 
e such that the turning points are given by 



pM 
pM 



(15) 
(16) 



where < e < 1. These two constants describe the 
geometry of the orbit, just as in Keplerian orbits: p is a 
measure of the radial extension of the orbit, while e is a 
measure of its deviation from circularity. These constants 
can be related to E and L by letting f = in Eq. (fT3|) . 
which leads to 



E^ 



(p-2-2e)(p-2 + 2e) 



p{p' 
p'^M^ 
p — i — e' 



3-e2) 



(17) 
(18) 



Continuing to exploit the analogy with Keplerian or- 
bits, we introduce a parameter x that runs from to 27r 
over one radial cycle, such that r{x) takes the elliptical 
form 



pM 



1 + e cos(x — w) ' 



(19) 



where w is the value of x periapsis, referred to as 
the argument of periapsis. The radial component of the 
velocity is hence 



r'ix) 



pMe sin(x — w) 
[l + ecos(x — w)\ 



(20) 



where a prime henceforth indicates a derivative with re- 
spect to X- 



From these results we can relate the parameter x to 
sing ^ — which yi' 

p3/2M(p-3-e2)i/2 



the proper time r using ^ — which yields 



dr _ 

d-X (p — 6 — 2e costi)i/2(i 4- e cosu)2 ' 
where we have introduced the variable 
V = X — w 



(21) 



(22) 



for brevity. Along with Eqs. ((T2]), HT]), and p^ . 

this leads to the following parametrizations for t{x) and 



^ix) 
<P'ix) 
tix) 
fix) 



= $ 



(l>'{x)dx, 



P 



p — 6 — 2e cos V ' 

T+ rt'{x)dx, 



p^M 



(p — 2 — 2ecosti)(l H- ecosti)^ 

^ / (p-2-2e)(p-2H-2e) 
y p — 6 — 2e cos v 



(23) 
(24) 
(25) 

(26) 



where we have defined the constants T and $ as the 
values of t and (j) at periapsis, respectively. 

Our parametrization of bound geodesies consists of 
Eqs. nil), (Uni), and (US!)-®- We see that a geodesic 
is uniquely specified by the orbital elements — 
{p, e, ly, T, $}. The principal elements p and e deter- 
mine the spatial shape of the orbit and are equivalent 
to specifications of energy and angular momentum; they 
determine the choice of geodesic. The positional ele- 
ments w, T, and $ determine the spatial orientation 
and time-dependence of the orbit; they determine the 
starting point of the particle on the selected geodesic. 
All together, the specification of the orbital elements is 
equivalent to the specification of initial values for the po- 
sition and velocity of the particle. We need three initial 
positions for a planar orbit, and we need two initial ve- 
locities (three minus one, by virtue of the normalization 
condition on the velocity vector); this counting matches 
the number of orbital elements. 

We note that our choice of orbital elements is closely 
related to Mino's in Ref. [l3|. When the orbital motion 
is restricted to the equatorial plane of a Kerr black hole, 
Mino uses the principal elements E and L and positional 
elements that are identical to our w, T, and $. To use 
(p, e) instead of {E, L) is mostly a matter of taste; we 
believe that the set (p, e) is more useful than {E, L) be- 
cause it gives a simpler parametrization, and because p 
and e are geometrically more informative. In the follow- 
ing subsection we will deviate more strongly from Mino's 
parametrization: for reasons that will be explained, we 
shall avoid directly evolving the elements T and <i>. 
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All the equations presented in this section remain valid 
for a perturbed orbit, with the exception of Eqs. (jlSp and 
(fT6| . which lose their meaning. The alteration that we 
shall make to account for the perturbation is that in each 
equation, the orbital elements will become functions of x- 



C. Evolution equations 

If we restrict the perturbing force to lie in the plane of 
the orbit, and assume that the orbit remains bound, then 
the geodesies described in the last section form a suffi- 
cient set of osculating orbits. Using our parametrization 
of these geodesies, along with the results of our general 
analysis in Sec. Ill Al we can now find evolution equa- 
tions for the orbital elements. Multiplying both sides of 
Eq. (Uni) by ^, we find 







dr 1 dr 






0, 


(27) 




dr , ^ 
dp^ 


de dw 


w' 




dt , 


dt , 


dt , 








(28) 


dp^ 




+ IT^' + 
dw 


r 




0, 


d± , 
dp 


de 


d(t> , 
dw 


$' 




0. 


(29) 


rom Eq. (fTTj) 


we find 










di 


, di 


, di , 








(30) 


dp^ 


de 


aw 




fr' 




dr 


, dr 






fr 




(31) 


dp^ 


^d-e' 


dw 




f 

■> 


d± , 
dp 


. d^ 
de 


dw 






1 


(32) 



The orthogonality condition /"ia = allows us to 
remove one component of Eq. ([TT]) from the set of equa- 
tions; we use this freedom to remove Eq. ([30)) . The re- 
maining equations decouple into a closed system of ordi- 
nary differential equations for p, e, and w and two auxil- 
iary equations for T and $. We shall find that the evolu- 
tion equations for p, e, and w are simple. The equations 
for T and <i>, however, are not: Factors such as ^ in 
Eqs. ((28)) and ((29)) introduce elliptic integrals of the form 

J w dp 



integrals would have to be evaluated at each time-step in 
a numerical evolution, and they would create an exces- 
sive computational cost. Additionally, the integrals gen- 
erally grow linearly with Xi and this produces terms in 
r(x) and <&(x) that grow quadratically with x, as well as 
terms that oscillate with a linearly increasing amplitude. 
Such terms greatly confuse both numerical and analyti- 
cal descriptions, and they are largely an artefact of our 
parametrization. (This statement applies also to Mino's 
parametrization [13] ■) We note that similar (though less 
severe) difficulties arise also in the method of osculating 
orbits in Newtonian celestial mechanics; refer for example 
to the discussion on pp. 248-250 in the text by Beutler 
[lol |. In the Newtonian context, alternative orbital ele- 
ments are typically selected so as to overcome these prob- 
lems. With no obvious choice of alternative elements in 
the relativistic context, we opt instead to directly evolve 
the coordinates t and (f> rather than the elements T and 

Our phase space thus consists of {p,e,w,t,(j)}. This 
choice of phase space does not allow an easy separation 
of perturbative from geodesic effects in the evolutions of t 
and (j), nor does it allow a clean separation of conservative 
from dissipative effects. But it is overwhelmingly more 
convenient than the alternative choice {p, e, w, T, $}. If T 
and $ are required in an application, they may be found 
as, e.g., T = t — J^t'{x)dx- This may be necessary if 
initial conditions are required on an osculating orbit, or 
if one wishes to fully isolate perturbative effects. 

Solving for w' from Eq. ((27)) . and noting that = 
— r', we find 



\ f dr 



I ^ dr ^, 
dp^ de 



(33) 



Substituting this into Eqs. ((30)) and ((32)) . we can solve 
for p' and e' to find 



Ce{c^)Cp{r) ^ Ce{r)Cp{cj)f ' 
Cp{r)f'i> - Lp{4>)r , 



(34) 
(35) 



Iw ^(x)^X iiito the expressions for T' and These 



Ce{(t))Cp{r) ^ Ce{r)Cp{(j)) ' 
where £a(a^) = ff + Explicitly, the results are 
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P = 



2p'^/^M'^{p - 3 - e2)(p - 6 - 2e cost;) (p _ 3 _ cos^ v) ^ 2p^Me{p - 3 - e^) sin w 

(39 - 6 + 2e) (p - 6 - 2e) (1 + e cos v)* ^ (p - 6 + 2e)(p - 6 - 2e)(l + ecos?;)^ 

p5/2M2 (p - 3 - e^) { (p - 6 - 2e2) [(p - 6 - 2e cos u)e cos t; + 2(p - 3)] cos v + e(p2 - lOp + 12 + 46^) } 
(p - 6 + 2e) (p - 6 - 2e) (p - 6 - 2e cos t;) 1/2(1 + e cos u)"* 
p^Af (p - 3 - e2)(p - 6 - 2e2) sin w 



(p - 6 + 2e)(p - 6 - 2e)(l + ecosw)2 
p5/2j\,f2(p _ 3 _ e^) {(p - 6) [(p - 6 - 2ecosw)ecosw + 2(p - 3)] - 4e^ cosv} sinw 
e(p — 6 + 2e)(p — 6 — 2e)(p — 6 — 2ecosi;)i/2(l + e cosw)** 



p2M(p - 3 ~ e2) [(p - 6) cos v + 2e] 
e(p — 6 + 2e)(p — 6 — 2e)(l + ecosu)' 



r,(36) 

(37) 
(38) 



These equations could be rewritten in any number of 
ways, in terms of alternative linear combinations of /*, 
/'', and /"^j by using the orthogonality relation /ai" ~ 0, 
which has the explicit form 



0. 



(39) 



The result of such a rearrangement might in fact be sim- 
pler, but it may also be ill-behaved from a numerical 
point of view. One such alternative combination is given 
in Appendix [Bl 

Our first formulation of the method of osculating or- 
bits is complete. We have first-order evolution equa- 
tions for each one of the dynamical variables in the set 
{p, e, w, t, (j)}] the equations for t and (f> were obtained in 
the preceding subsection, and for convenience they are 
reproduced here: 



t' = 



p2M 



(p — 2 — 2ecosti)(l -t- ecosv)2 

^ / (p-2-2e)(p-2-f-2e) 
V n — 6 — 2e cos V 



P 



p — 6 — 2e cos V 



(40) 
(41) 



Equations ([571) . ^^'^ ([55]) form a complete set of 

equations for p(x), e(x), and w(x); once these functions 
are known, t{x) and (j){x) can be obtained from the re- 
maining two equations. We recall that w = x — w{x)- 

One may note that w' diverges as e ^ 0. This corre- 
sponds to the fact that w loses its geometric meaning for 
circular orbits. To overcome this difficulty we can again 
follow celestial mechanics and define alternative orbital 
elements a = esinw and /3 = ecosw. The radial coordi- 
nate in terms of these elements is 



pM 



1 + + 



(42) 



where \1/ = a sin x and Q = (3 cos x are introduced 
for the sake of brevity in later expressions. While a 
and (3 do not possess a clear geometric meaning, which 
limits their usefulness for generic orbits, they do allow 
one to analyze small-eccentricity or quasi-circular orbits. 
Their evolution equations can be easily calculated as 
a' — e' sin w + ew' cos w and /?' = e' cos w — ew' sin w. Us- 
ing the identities e cos v = a sin x + /3 cos x and e sin v — 
/3sinx — acosx to simplify the results, we find 
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p5/2M2(p-3-a2 



a/3 cos 2x + - (q^ - ) sin 2x 



Vp~~6~2(¥TTI)((p - 6)2 - 4(a2 + fj2)^^{i + m + 
+ [2{p - 3) + (p - 6)(1' + f7) - 2{-^ + nf] [{p - 6) cosx - 2/3(«' + f7)] + /3 [p^ - lOp + 12 + 4(a2 + /32)]j 
p'^M{p - 3 - - /32) [(p - 6 - 2/3^) sinx + 2a(l + f^)] J'' 



((p - 6)2 - 4(a2 + /32))(1 + * + f])2 



(43) 



95/2m2(p-3 



/32)/^ 



X a/3 



a/3 cos 2x + - (a^ - /5^ ) sin 2% 



^39-6-2(* + r2)((p - 6)2 - 4(a2 + /32))(1 + * + fi)^ 

+ [2{p - 3) + (p - 6)('J' + f7) - 2{-^ + nf] [{p - 6) sinx - 2a{^ + n)] + a [p^ - iQp + 12 + 4(a2+/32)]J 

p'^M{p-3-a^ - (3^) [(p-6-2a2)cosx + 2/3(l + *)] f 
((p - 6)2 - 4(a2 + /52))(1 + ^ + n)^ ■ 

To evolve our full system we must also express p', t' , and </>' in terms of a and /3: 



(44) 



2p^/2AfVP - 6 - 2(«' + f^)(p - 3 - a2 - /32)(p - 3 - («- + 17)^ 
[(p - 6)2 - 4(a2 + /32)](1 + * + r!)4 

2p^M(p- 3- a^ - /32)(/3sinx- acosx)/"" 
[(p - 6)2 - 4(a2 + /32)](1 + * + f7)2 ' 

p2 M v/(p-2)2-4(a2+/32) 
(p - 2 - 2(* + r^)) Vp~~6~2(¥TT^(1 + * + S1)2 ' 



p- 6 -2(1' + 17)' 



(45) 
(46) 
(47) 



r 



This is our second formulation of the method of osculat- 
ing orbits. The first formulation involves shorter equa- 
tions, but it becomes ill-behaved when e is small. The 
second formulation is well behaved, but it involves longer 
equations. 



with the 2.5PN equations of motion for each one of the 
two bodies. Within the center-of-mass frame the relative 
motion of the two bodies is governed by the closed system 
of equations [13] 



III. POST-NEWTONIAN BINARIES 



A. Hybrid equations of motion 



dt 



-IT . (48) 



We now move on to a concrete application of our 
method by considering the post-Newtonian binary sys- 
tem introduced in Sec. II Bl This system consists of two 
gravitationally-bound bodies of mass mi and m2, with 
equations of motion derived to 2.5PN order in a post- 
Newtonian expansion; because we are interested in self- 
force effects, we take the ratio mi/m2 to be small, and we 
neglect the spin of the bodies. In this section we explain 
how such a system can be analyzed with our method of 
osculating orbits. 

Our analysis is based upon the hybrid equations of 
motion presented in Ref. T^j- These equations begin 



where xf^ = — is a Cartesian spatial vector from 
TO2 to mi in harmonic coordinates, = Sabx"'x^ is the 
square of the vector's Euclidean magnitude, t is a har- 
monic time coordinate, and M = mi + m2 is the total 
mass of the system. The functions A and B depend only 
on the total mass A/, the reduced mass /i — TO1TO2/M, 
and the relative coordinates and velocities. They can be 
written as A — Am + ^A and B = Bm + eS, where 
e = [ijM and terms with a subscript M are independent 
of /I. The /x-dependent terms are quadratic in e, and 
they can be further decomposed into post-Newtonian or- 
ders as A = Jii + ^2 + ^2.5 and B = Bi+ B2 + -82.5- 
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Explicitly, these have the form 
, M 



Am - 1-4— + 1.2 + 9 ( — ) -2— i:^) , (49) 
Th \rhj rh \dt ' 



eAi = -e 



2. 6v + - — — 

Th 2\ dt 



(50) 



M 



-ri — ] +(3-4e)t;4--(13-4e)— z; 



drh 



15, 



(25 + 2e) 



dt J ' 8 

2 - 



rh \ dt 



2.5 
Bm 

eB2 



8 Mdrft, 
5 r/,. 



17 M 

3 Th^ 



, ,4-2— 
dt V Tft 



2e 



eB 



2.5 



dvh 
t 

drh 
dt ' 
1 c;rft_ 
2^ dt 



8 M 
5 r,i 



(15 + 4e)u2 - 

-3(3 + 2e) 

+ 3— 

rh 



M 

(41 + 8e) — 

rh 



drh 
dt 



(51) 

(52) 

(53) 
(54) 

(55) 
(56) 



where = ^ab-^-^t square of the velocity vector 

in harmonic coordinates. 

The hybrid equations are inspired by the fact that 
when e = 0, Eq. becomes identical to a 2PN expan- 
sion of the geodesic equation in a Schwarzschild space- 
time with mass parameter M. Building on this fact, 
Kidder, Will, and Wiseman [l3| replaced Am and Bm 
with their exact geodesic expressions As and Bs in the 
fictitious Schwarzschild spacetime. In other words, the 
hybrid equations of motion are given by Eq. (I48p after 
substituting A = As + eA and B = Bs + eB, where 



A, 



Bs 



1 - M/rh 
(1 + M/rhf 

2- M/rh M 
^ 1 - MVr2 'V^ 

4 - 2M/rh drh 
1-A'P/rl dt 



drh 
dt 



(57) 
(58) 



The resulting equations are accurate to 2.5PN order, but 
in the test-mass limit mi ^ they exactly describe 
the orbit of the test mass in the Schwarzschild space- 
time of the other body. These equations form an ideal 
test case for our method of osculating orbits because, 
besides their relative simplicity, they explicitly split into 



geodesic terms and perturbation terms. This allows us to 
construct osculating orbits as geodesies in the fictitious 
Schwarzschild spacetime of mass M. We can then easily 
derive the perturbing force from the terms A and B. 

The first step in this process is to write the equations 
of motion in plane polar coordinates {rh,(j)), which are 
defined by x\ — rhCos{4>) and x| = r/isin((/)). In terms 
of these coordinates, Eq. becomes 



frh 
dt^ 

dt^ 



M 

M d(t> 
~tB—— 

rl dt 



A + B- 



.drh 



dt 

2 drfi d(j) 
rh 



rh 



dt 



dt dt ' 



(59) 



(60) 



The harmonic coordinates used here are related to 
Schwarzschild coordinates by the simple transformation 
r/j = r — M . Since M is constant, the subscript h can 
be safely dropped within derivatives. Expressing r/j in 
terms of r, the above equations are transformed into 
Schwarzschild coordinates. 

We derive /" from these equations as follows. From 
Eq. ^ we have 



r 



dt^ 



dz^ dz'^ 



dz' 

n 



(61) 



Although we could calculate K{t) directly from its defini- 
tion, the result would be unwieldy. We instead use the 
equation of motion for t, 



dH 

df^ 

to replace k with 



dt 



J dz^ dz'^ J ._2 



Substituting this expression for k into Eq. 



dz' 



(62) 

(63) 
1]), we find 

(64) 



where 



/pa rf^"pt \dZ-dZ' 



dzf dz'^ 



dt^ 



(65) 



The subscript p refers to the fact that involves only 
the perturbative terms in d^z^/dt^. Indeed, a simple 
calculation based on the preceding equations for d^r/dt^ 
and d^(f)/dt^, as well as the Christoffel symbols obtained 
from the Schwarzschild metric, reveals that 



~dr 
eB — 
dt 



(66) 
(67) 
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Equation ((64)) determines and in terms of /*. 
The orthogonality condition ([5^ then aUows us to find 
all three components of the force. The result is 



P dt 



"p' ^ dt 



F2 



"p ' dt dt 



F 



rdrd± I 0/^p2 _ ( dr\'^\ 
"p dt dt ^ '^p\^ \dt) ) 



F2 - 



' dr_ 

V dt 



^"-Fr^{tf 



(68) 



(69) 



(70) 



Substituting into the above results, and using the nor- 
malization condition —1 = = — Ft^ + F^^r^+r^f/)^, 
leads to 



F 



eMt-* d<t) 
rl 'dt 



^F-^{dr/dt)A 



FB 



{dr/dt)BY 
}• (71) 



Since /* is not required in our formalism, we will not 
provide an explicit expression for it. 

We can recast these equations in a form analogous to 
that of Eqs. ((591) and ((60)) . 



A + B 



,dr 
di 



f 



'r^ dt' 



by defining A and B as 



A = 



B 



(1 - M/rf 
(1 - M/rf 



1 dr ~ 

A 

F dt 



FB 



(72) 
(73) 

(74) 
(75) 



The factors of t convert the "time" variable in the accel- 
eration from coordinate time to proper time; this is given 

by 



for A and B into Eqs. ((74)) and ([75]); the relevant equa- 
tions are listed near the beginning of Sec. IIII Al In these 
equations we must make the substitution rh = r — M, 
and convert t-derivatives into ^-derivatives by employ- 
ing Eq. ((26)) . In these final forms, the expressions for 
and are ready to be inserted within the evolution 
equations for the orbital elements. 



B. Results 

1. Adiabatic, secular, and radiative approximations 

We are primarily interested in determining the types 
of errors introduced by the adiabatic and radiative ap- 
proximations. We should first clarify the meaning of 
these approximations. The basis of both approximations 
in the context of osculating orbits is the separation of 
orbital elements into secular and oscillating parts, i.e. 

— + /^(. . The particular adiabatic approximation 
that we are concerned with, which we have titled "secular 
approximation," is one which eliminates the oscillations 
and keeps only the secular behavior; that is, it uses an 
approximate orbital evolution with J^j^ — I^^. A ra- 
diative approximation uses only dissipative terms in the 
perturbing force, with orbital elements I^, with the hope 
that the secular part ifg^^. of this evolution reproduces 

sec 

Unfortunately, these general definitions are somewhat 
ambiguous. We examine first the case of the secular ap- 
proximation. The main source of ambiguity associated 
with the general idea of removing oscillations is that it is 
not clear which oscillations are intended to be removed. 
For example, in the formalism presented in this paper, 
removing the oscillations with respect to x will not re- 
move the oscillations with respect to t, and vice versa. 
This failure is caused by the zeroth-order (i.e. geodesic) 
oscillations in time as a function of x- Consequently, a 
secular evolution defined by an average over the orbital 
parameter x, such as 

lL = {I^)x^^ r^I^{x')dx', (77) 



t' 



F - F-^{dr/dtY - r^{d(j)/dtf 



(76) 



where, we recall, F = 1 — 2M/r. The factors of 
1/(1 — Af/r)2, on the other hand, convert from harmonic 
coordinates to Schwarzschild coordinates. One could in- 
corporate these factors into each Ai and Bi and then 
re-expand these in powers of M /r to find new expres- 
sions for Ai and Bi, neglecting terms of 3PN order and 
higher; but since the hybrid equations already introduce 
errors above 2.5PN order, doing so is unnecessary. Thus, 
for simplicity we shall use the force in its above form. 

The final expression for the perturbing force is ob- 
tained by substituting the post-Newtonian expansions 



will differ from that defined by an average over time, such 
as 



rx+-jAdi^, 

J-y-n dy 



(78) 



A precise definition of a secular approximation would 
have to specify which averaging procedure is to be se- 
lected. 

A second source of ambiguity concerns the choice of 
initial conditions. We desire that our secular evolution 
reproduce the average of the true evolution, and this 
means that in general, the initial conditions placed on 
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FIG. 1: Comparisons of true orbits (solid black curves) and ra- 
diative approximation orbits (dashed green curve) with iden- 
tical initial conditions and with a mass ratio fi/M = 0.01. 
In each case the two orbits begin at periapsis and are termi- 
nated at the same final time. Upper plot: highly eccentric 
orbits with po = 50 and eo — 0.9. At the end of the sim- 
ulation the approximate orbit lags behind the true orbit by 
approximately one-half radial cycle out of a total of fifteen. 
Lower plot: quasi-circular orbits with identical initial condi- 
tions Po = 10 and eo = 0. Again, the approximate orbit lags 
behind the true orbit. 



the approximate solution will have to differ from the ex- 
act initial conditions. This is because the exact solu- 
tion contains the secular approximation plus oscillations, 
and the oscillations may not vanish at the initial time. 
Identifying the correct initial conditions for the approx- 
imate evolution therefore requires knowledge of the os- 
cillations; in the absence of such information — that is, 
when the exact solution is not known — the initial condi- 
tions remain unknown and the procedure is ambiguous. 
The ambiguity persists even when the exact solution is 
known, because it is then inherited from the first source 
of ambiguity, the question as to which oscillations are to 
be removed. The ambiguity associated with the initial 
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FIG. 2: The same eccentric orbit as shown in Fig. [T] but 
now using time-averaged initial conditions for the radiative 
approximation. In this case the approximate orbit is indis- 
tinguishable from the true orbit on the timescale of the plot 
(fifteen orbital cycles). 

conditions is lifted only when the averaging procedure 
is selected, and when the exact solution is known; the 
approximate initial conditions are then calculated by av- 
eraging the exact evolution over the first radial cycle. For 
example, in the case of an averaging over x we would set 

4l(0) = (27r)-i/f /^(X)cix. 

We shall not pursue a detailed exploration of the ambi- 
guities associated with the secular approximation in this 
paper, although they are quite important; they are the 
focus of a companion paper [l7j . Our focus here will be 
instead on the limitations and ambiguities of the radia- 
tive approximation. As was indicated previously, a radia- 
tive evolution switches off all conservative terms in the 
perturbing force (Ai = = -Bi = i?2 = 0), and retains 
only the radiative terms at 2.5PN order (A2.5 7^ and 
^2.5 7^ 0). This approximation is logically distinct from 
adiabatic approximations in general, but the hope formu- 
lated in the literature (for example in Refs. J.8, 19, 20]) 
is that the radiative evolution will reproduce the secular 
changes of the orbital elements. We shall see that while 
the radiative approximation captures the secular changes 
in p(x) and e(x), it fails to account for secular changes 
in w(x), i(x)i and (/>(x); this conclusion confirms and 
extends those of our previous work 

In addition, the radiative approximation is subject to 
the same ambiguities regarding the choice of initial con- 
ditions as the secular approximation. Writing the radia- 
tive evolution as the sum of its secular and oscillatory 
parts, l^{x) = //'ggc + ^r^osc I shall consider three pos- 
sible candidates for J/'(0). The first is //^(O) = /^(O), 
the exact initial data that is selected for the true evo- 
lution of the orbital elements under the action of the 
full perturbing force. The second is 1^^^^ — (-^'^)x(0)i 
the x-averaged initial data, which identifies the initial 
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FIG. 3: The principal element p and positional element w as 
functions of time for a complete inspiral, beginning with the 
initial conditions of the eccentric orbit in Fig. [1] In each plot 
the true curve is in solid black, the radiative curve with the 
same initial conditions is long-dashed in green (the uppermost 
curve in the p plot) , the radiative curve with x-averaged initial 
conditions is short-dashed in blue (middle curve in p plot), 
and the radiative curve with time-averaged initial conditions 
is dotted in black (lowest curve in p plot). The insets display 
the early behavior of the curves, covering the same range of 
time as in Fig. [1] 



secular part of the radiative evolution with the initial 
X-averaged part of the true evolution. The third choice 
is /i^sec(O) = (-^'^)t(O)' the t-averaged initial data, which 




10000 20000 
X 



30000 



FIG. 4: The principal element p as a function of the orbital 
parameter x- The curves are as described in Fig. [S] The 
radiative curves do deviate secularly from the true curve, but 
the errors are too small too appear on the scale of the graph. 



identifies the initial secular part of the radiative evolution 
with the initial t-averaged part of the true evolution. (We 
note that for both the second and third choices of initial 
conditions, the initial value /i^(0) is not fixed by J^sec(O) 
alone, since we also require the initial value of if^sc- 
though we do not have a priori access to this oscillatory 
part, we can assign it an approximate initial value based 
on the results of the radiative evolution with exact initial 
conditions. This introduces a negligible error, since the 
oscillations in the radiative evolution are extremely small 
in practice.) These three choices of initial data are dis- 
tinct, and they lead to different evolutions. We shall see 
that the accuracy of the evolution (relative to the true 
evolution) depends strongly on the choice of initial data. 



2. Orbital evolution 

A typical inspiral of interest for LISA will form in a 
highly eccentric state. Over the course of the inspiral the 
system will emit gravitational radiation carrying away 
energy and angular momentum, shrinking and circular- 
izing the orbit over time. Thus, the inspiral will evolve 
from a highly eccentric orbit to a quasi-circular one, and 
it will end in a rapid plunge. We shall now determine 
the validity of the radiative approximation for this class 
of orbits. Since our perturbing force is valid only in the 
post-Newtonian regime, we always ensure that ^0.1. 

The general limitations of the radiative approxima- 
tion are demonstrated in Fig. [1] which displays the spa- 
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tial trajectories of a highly eccentric orbit and a quasi- 
circular orbit, along with corresponding radiative approx- 
imations. In each case the true and approximate orbits 
are terminated at identical final times, at which point 
the radiative approximation lags behind the true orbit. 
With a mass ratio of /i/Af — 0.01, this dephasing of the 
two orbits is noticeable after only fifteen radial cycles 
in the eccentric case, while several dozen revolutions are 
required in the quasi-circular case. Since the dephasing 
is apparent before any non-geodesic precession occurs, 
we interpret its cause to be conservative effects in the 
time-dependence of the orbit. That is, the error in 
dominates over the errors in u'(x) and ^(x)j such that 
the particle lies at the wrong spatial point at a given 
time, even before and 0(x) have deviated signifi- 
cantly from the true orbit. 

For the plots in Fig. [1] we have chosen exact initial 
conditions l^ifS) for the approximate orbit. By choosing 
averaged initial conditions we obtain better results in the 
eccentric case: as shown in Fig. [2l using time-averaged 
initial conditions (/"^)t(0) eliminates the dephasing on 
the timescale of the plot. Using x-averaged initial con- 
ditions (/'^)x(O) results in a smaller improvement, as we 
will discuss below. However, in the quasi-circular case all 
initial conditions fare equally well. 

The evolution of the orbital elements over a complete 
inspiral, beginning with the initial conditions of the ec- 
centric orbit in Fig. [T] and continuing to quasi-circularity, 
is displayed in Fig. [31 Insets in the plots display the 
same range of time covered by Fig. [T] The orbit stops 
before the final plunge of the small body into the large 
black hole. There are two reasons for this truncation. 
First, our method of osculating orbits cannot cover the 
final plunge, because of the underlying restriction that 
the orbit must be bounded between a minimum radius 
pM / [1 + e) and a maximum radius pM/{l — e); this is 
reflected mathematically by the condition p > 6 + 2e, 
which is violated during plunge. Second, we should in 
any case leave this portion of the orbit alone, because 
the velocities and fields therein are highly relativistic; 
in this regime the post-Newtonian expansion of the per- 
turbing force becomes inaccurate. In Fig. [3] we display 
results of the numerical evolution for the principal ele- 
ment p and positional element w only; the evolution of 
e is qualitatively similar to that of p. It is worth noting, 
however, that the eccentricity never quite reaches e w 0; 
instead, quasi-circularity is manifested by the condition 
X — w ~ 0, which equally well ensures that r' ~ 0. This 
observation agrees with the results of Ref. f23| . 

The results for all three choices of initial conditions 
are plotted in Fig. [3] As we see from these plots, the 
radiative approximation qualitatively matches the true 
secular evolution for the principal element p, but neglects 
all secular changes in the positional element w. This is 
the expected result. However, we also see that the ra- 
diative approximation deviates from the true evolution 
even for the principal element. The extent of this devia- 
tion depends on the choice of initial conditions, with the 
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FIG. 5: The orbital phase 0, with curves as described in 
Fig. |3l The scale of the plot suggests that the radiative evo- 
lution with time-averaged initial data (the uppermost curve 
in dotted black) gives an accurate approximation of the true 
evolution. The vertical scale, however, is large, and this is 
a false impression. At the late time t/M — 2.345 x 10^, the 
error in phase is = 4520 rad for the exact initial data, 
A<j!> — 1830 rad for the x-averaged initial data, and Ai^ = 655 
rad for the f-averaged initial data. This last choice fares best, 
but its accuracy is poor over a complete inspiral. 



time-averaged initial conditions faring the best and exact 
initial conditions the worst. 

An essential aspect of our results is that the errors 
in the principal elements produced by the radiative ap- 
proximation are mostly due to errors in t{x)- As we see 
in Fig. 21 the errors almost completely vanish when the 
principal elements are plotted as functions of x! signif- 
icant errors arise only in the conversion between x ^^nd 
t. 



3. Errors in orbital phase 

The errors in which we are most interested are errors 
in orbital phase, since they will lead directly to errors 
in the phase of the emitted gravitational radiation. Fig- 
ure [5] displays the phase 4> versus time, again using all 
three choices of initial conditions for the radiative ap- 
proximation. Once again we see that the time-averaged 
conditions produce the smallest error, for the same rea- 
sons described in the previous section. 

Figure [6] shows the dependence of the dephasing Acj) = 
^ '/'rad on the parameters of the problem. We plot 
the dephasing for a "radiation-reaction" time defined 
by P — 0.9poi rather than a complete inspiral, since 
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FIG. 6: The difFerence in orbital phase (p between the true 
orbit and approximate orbits after a radiation-reaction time 
(defined by p — > 0.9po on the true orbit). Open squares indi- 
cate resuhs for identical initial conditions, open triangles for 
matching ^-averaged initial conditions, and open circles for 
matching f-averaged initial conditions. Top: dephasing as a 
function of the mass ratio n/M, with fixed true initial values 
po = 50 and eo — 0.9. The dephasing becomes /^-independent 
for sufficiently small fi, when second-order effects become neg- 
ligible. Middle: dephasing as a function of initial value po, 
with fixed eo — 0.9 and i-i/M = 0.1. Bottom: dephasing as 
a function of initial eccentricity eo, with fixed po = 50 and 
jj,/M = 0.1. The error in the case of time-averaged initial 
conditions is approximately independent of e. 



gravitational-wave data analysis may require only a por- 
tion of a complete inspiral. We see that the dephasing is 



independent of for sufficiently small values of fj,. This 
is an expected result, since the radiation-reaction time 
at leading order in fi varies as l//i, while the rate of de- 
phasing varies as /i, leading to a net cancellation in the 
total dephasing. However, terms in the perturbing force 
that are quadratic in fi alter this result when /i/M is suf- 
ficiently large. Somewhat surprisingly, these quadratic 
terms actually serve to decrease the dephasing, lowering 
the impact of conservative terms in the force. 

As expected, the dephasing decreases at lower values of 
e, although the eccentricity seems to have negligible im- 
pact in the case of time- averaged initial conditions. Also 
as expected, the dephasing varies as p^^^, regardless of 
initial conditions. This scaling follows from the form of 
the post-Newtonian force: the leading-order conservative 
term enters at IPN order, which scales as a correc- 
tion to Newtonian gravitation, while the leading-order 
dissipative term enters at 2.5PN order, which scales as 
a correction. The dephasing is governed by the 

relative strength of the conservative terms, leading to a 
scaling of p~^/p~^^'^ — p^^^. 

In all cases the time-averaged initial conditions yield 
the best results. Indeed, the efficacy of these initial con- 
ditions is almost surprising. One way of understanding 
their impact is to examine the insets in Fig. [31 Peaks in 
the true curve correspond to the short periods of time 
near periapsis, while relatively flat regions correspond to 
the long periods of time around apoapsis. Thus, choos- 
ing exact initial conditions matches the true and approx- 
imate orbits for the minimal amount of time, as well as in 
the region of strongest fields, leading to the largest possi- 
ble deviation. Choosing time-averaged initial conditions 
matches the orbits near apoapsis, for the longest time and 
with the weakest fields, leading to the least possible devi- 
ation. The x-averaged initial conditions are then in some 
sense the average of all the incorrect choices. An implica- 
tion of this is that in some circumstances the x-averaged 
initial conditions could turn out to be even worse than 
the exact initial conditions. For example, choosing exact 
initial conditions at apoapsis would closely approximate 
the time-averaged initial conditions, which would then 
fare much better than the x-averaged initial conditions. 

Wc can explain the long-term impact of the initial 
conditions by considering the time-dependence of an or- 
bit. The secular time function (t)(x) can be written 
in terms of the orbital period P{x) = J^^^t'{x)dx as 

= Jo Pix)dx- As we see from the insets in Fig. [3l 
the changes in initial conditions bring the initial orbital 
period of the radiative approximation closer to that of the 
true orbit; and as we would intuitively expect, the time- 
averaged initial conditions best reproduce the initial tem- 
poral period. This correction, SP, to the initial period 
then induces a long-term correction to (t) (x) of the form 

X • SP- (Such an effect can be calculated expfi citly 
for the electromagnetic self- force considered in Ref. [rn |: 
refer to Eqs. (4.8) and (4.21) therein.) In essence, the 
time-averaged initial conditions carry information about 
the initial conservative correction to the true orbital pe- 
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riod, and they thus implicitly insert a conservative cor- 
rection into the radiative approximation. This serves to 
remind us that we would have difficulty choosing suit- 
able initial conditions for the radiative approximation if 
we did not have prior access to the true evolution. 

Regardless of the choice of initial conditions, one 
should note that the errors accumulated over a complete 
inspiral are much larger than those shown in Fig.[Sl (Re- 
fer to the caption of Fig. [5] for actual values.) Also, the 
plots of versus p and e are for e = 0.1, leading to 
a smaller dephasing than would occur if e were in the 
region of linear dominance. Thus, even if ideal initial 
conditions could be found without reference to the exact 
solution, the radiative approximation would generically 
fail over a complete inspiral. 



4- Gauge dependence 

As is well known, the gravitational self-force is a gauge- 
dependent quantity: it is not invariant under a change of 
coordinates ^ + S,^ , where is a "small" vector 
field. The equations of motion that we have used in this 
paper were calculated within the harmonic gauge of post- 
Newtonian theory, and the magnitudes of the conserva- 
tive effects that we have displayed refer to this particular 
gauge choice; different gauges would necessarily lead to 
different results. Indeed, Mino has argued in favor of 
constructing a physically meaningful "radiation-reaction 
gauge" in which the conservative effects of the self-force 
are set to zero over a finite radiation-reaction time, mak- 
ing the radiative approximation exact over that interval 
[13, nil ■ Mino has also argued that this gauge choice in- 
duces a change in initial conditions that partially absorbs 
conservative effects [31, and this statement agrees with 
our result that long-term conservative effects can be mim- 
icked by a small change in initial conditions. We would 
like to point out, however, that a rigorous construction 
and implementation of such a gauge choice have yet to 
be performed, and that the impact of making this choice 
on quantities other than the self-force has yet to be de- 
termined. 

It is known, for example, that in the harmonic gauge 
of post-Newtonian theory, the equations of motion con- 
tain both radiative and conservative terms, and that the 
gravitational potentials are well-behaved everywhere, ex- 
cept at the position of each (pointlike) body where they 
diverge with an expected power of m/r. What is the be- 
havior of the gravitational potentials in Mino's radiation- 
reaction gauge? The answer is not known, and it would 
be interesting to investigate the issue in post-Newtonian 
theory. For example, one could determine the effect 
on the potentials of making a coordinate transformation 
that would turn off some of the conservative terms in the 
equations of motion (those that depend on e in Sec. IIIA); 
would this spoil the behavior of the potentials near the 
bodies, or perhaps elsewhere in the spacetime? Such an 
analysis would be revealing, and it would give indication 



as to whether Mino's scheme is likely to be successfully 
implemented. 

We believe that the Lorenz gauge of the gravitational 
self-force problem, which is in close mathematical anal- 
ogy with the harmonic gauge of post-Newtonian theory, 
is also in close physical analogy: it produces conserva- 
tive terms in the self-force, and it produces gravitational 
potentials that are well behaved everywhere (except at 
the position of the orbiting body). Given the successes 
of post- Newtonian theory in its harmonic-gauge formula- 
tion, we feel confident that the Lorenz gauge is ultimately 
a better choice of gauge for the gravitational self-force 
problem, in spite of the presence of conservative terms 
in the equations of motion. We shall therefore defer our 
judgment on the advantages of Mino's radiation-reaction 
gauge, and reiterate the importance of the conservative 
terms in the harmonic-gauge (or Lorenz-gauge) self-force. 
Our conclusions, to be sure, apply within the confines of 
the post-Newtonian harmonic gauge. But we contend 
that our conclusions are in fact generic: Outside of a 
finely-tuned gauge choice, one should expect the con- 
servative part of the self-force to produce large secular 
effects. 



IV. CONCLUSION 

The first part of this paper was devoted to the devel- 
opment of a method of osculating orbits to integrate the 
equations of motion that govern bound, accelerated or- 
bits in Schwarzschild spacetime. The method involves the 
phase-space variables {p,e,w,t,(f)}, which are expressed 
as functions of an orbital parameter x] each variable sat- 
isfies a first-order differential equation, and knowledge 
of these variables is sufficient to determine the worldline 
in spacetime. Although the method is limited to situa- 
tions in which the force acts within the orbital plane, this 
limitation can be overcome; in addition, the force is not 
assumed to be small. We show in Appendix [X] that for 
large values of p, our equations reduce to the standard 
perturbation equations of Newtonian celestial mechanics. 
The method has many potential applications, including 
the important one of permitting an implementation of the 
gravitational self-force. Most immediately, it provides an 
attractive conceptual and mathematical foundation for a 
perturbative approach to weakly accelerated orbits. And 
furthermore, the method is easy to implement in practice 
in a numerical code. 

In the second part of the paper we applied the method 
of osculating orbits to the inspiral of a small body into 
a Schwarzschild black hole of much larger mass. The 
perturbing force was calculated on the basis of the hy- 
brid Schwarzschild/post-Newtonian equations of motion 
of Kidder, Will, and Wiseman [l3|, and its effect on 
the orbiting body was obtained by numerical integration 
of our evolution equations for the dynamical variables 
{p, e, w, (j)}. This approach is well suited to a study of 
the limitations and ambiguities of adiabatic and radiative 
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approximations, which was carried out next. Specifically, 
we have illustrated the importance of conservative effects 
in the time dependence of the orbit, and we have es- 
tablished the advantage of choosing time-averaged initial 
conditions for the approximated orbital elements. This 
problem differs in many respects from the fully relativis- 
tic self-force problem, but it nevertheless captures many 
of its essential features. Our conclusions, therefore, might 
be expected to hold in the fully relativistic case for most 
choices of gauge. 
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APPENDIX A: NEWTONIAN LIMIT 

Since our work extends the standard methods of New- 
tonian celestial mechanics, it is a worthwhile endeavor to 
show that our equations reduce to those for perturbed 
Kcplcrian orbits in Newtonian mechanics. In this Ap- 
pendix we derive the Newtonian limit of our expressions 
by expanding in powers p~^] since oc w^, this 

is equivalent to a post-Newtonian expansion. We shall 
first describe the general relationship between the New- 
tonian and relativistic perturbing forces. Next we shall 
show that our geodesic parametrization reduces to Ke- 
plerian ellipses and that our evolution equations for the 
orbital elements p, e, and w reduce to Gauss' perturba- 
tion equations of celestial mechanics. 

Substituting the Christoffel symbols of the 
Schwarzschild metric into the equations of motion ^ 
yields the following equations for the force: 



AI 



1 M , 

2" 



fcj) 



F4>\ 



1 2M , 



(Al) 
(A2) 
(A3) 



where F = 1 — 2M/r. The time-component of the force 
can be written in a more useful form using the orthogo- 
nality relation ([55)1 . 

These expressions for the relativistic force differ non- 
trivially from those in the Newtonian case. We define 
the Newtonian perturbing force per unit mass, via 
Newton's second law: 



(A4) 



where a; is a 3-vector representing the spatial coordinates 
of the particle and g = —^r is the Newtonian gravita- 
tional acceleration. For convenience we have defined the 
Newtonian acceleration as the second derivative of x with 



respect to proper time rather than coordinate time. We 
also define the radial and tangential components of the 
perturbing force via 



(A5) 



where r and cj) form an orthonormal basis in the orbital 
plane. Given these definitions, writing x in polar coordi- 
nates (r, 0) leads to 



pip — 



2fd 



(A6) 
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Comparing the Newtonian and relativistic expressions 
for the perturbing force, we see they are related by the 
equations 



f = ^"■ + r(l-F)02 
M 



p4> 



FP -f F-^r^ - I] 



(A8) 
(A9) 



Thus, differs from F"^ by relativistic corrections, while 
/■^ differs from F"^ only by a factor of the orbital radius. 

We next consider our parametrization of geodesies. 
From Eqs. ([Ml), JIS]), and one trivially finds the 
leading-order terms in cj)' , t' , and x to be 
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e cos(x — w)]'^ ' 
ecos(x — w)]^ 



P- 



Thus, in the Newtonian limit we have 
and the resulting parametrization 

pM 



e cos — w 



dt 



p^/^M 



(AlO) 
(All) 

(A12) 

X and t = T 

(A13) 
(A14) 



In terms of the orbital elements, we see that w — ^ 

in the Newtonian limit. This corresponds to the loss 

of one degree of freedom, as we would expect from the 

fact that t in Newtonian physics is a universal parameter 

rather than a coordinate. We can also easily find that 

the energy and angular momentum per unit mass reduce 
1 — ^ 

io E = 1 — ^-j^ and L — respectively. The first 

term in E is the rest energy of the particle, while the 
second term is the Newtonian energy — 

With the exception of the inclusion of the rest mass, 
the above results are standard Keplerian relationships. 
Thus, our equations for the orbital elements should re- 
duce to those for perturbed Keplerian orbits. Substi- 
tuting Eqs. (|M0)) - (|M2| into Eqs. ([SHD, jM]), and ((X9|) . 
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we find the leading-order expressions for the perturbing 
force: 



_ pr 

A. F'f' 

= 

f = 



r 

e sin( 



(A15) 
(A16) 

VP 



These results allow us to expand Eqs. ([36)1 . (I37)l . and 
([38]) to find the leading-order expressions for the orbital 
elements: 



dp 
lit 
de 
'dt 



2p 



3/2 



1 + ecos(0 — w) 

e + 2 cos((/i — w) + e cos^((^ — w) ~ 



(A18) 



dw 



1 + e cos(0 — w) 
+Vp sin((/) - w)F'' (A19) 
ypM^/^ sin((/) - w) [2 + e cos(0 - w)] 
e 1 + e cos(0 — w) 



cos(0 — w)F^. 



(A20) 



These are Gauss' well known perturbation equations. 



APPENDIX B: EVOLUTION EQUATIONS FROM 
KILLING VECTORS 



It is possible to derive Eqs. ([36l) -([38 l) for the deriva- 
tives of the osculating elements from Eq. (fTU]) and the 
Killing vectors of the Schwarzschild spacetime, without 
reference to Eq. (fTTjl . Although this derivation is equiv- 
alent to that given in Sec. Ill Ci its physical significance 
is more intuitive. We begin by defining energy and an- 
gular momentum (per unit mass) as E = — ^^-jia and 



L = £."^fa, where ^(t) = and ^(0) = ^ are Kilhng 
vectors corresponding to the spacetime's invariance un- 
der time translations and spatial rotations. From these 
definitions we find 



E 
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The first term on the second line vanishes due to the an- 
tisymmetry of S^a-p for any Killing vector i^, and the final 
line then follows from the equation of motion z'^z^-a — 
f^. An analogous result holds for L. From the definitions 
of ^(t) and 1^(0) we then find 



E Ff, 
L = r^f^. 



(B2) 
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These results can be used to find e and p using 
Eqs. (Ull) and ([TH]), which define E{p, e) and L{p, e). Us- 
ing these relationships, we write E = ^p + and 
L — ^p + §^e, which can be rearranged to find 
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dL dE 
de dp 



(B4) 
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The equation for w can then be found from Eq. (flOl 
which leads to Eq. or 



\ f dr , dr 



de dp^ 

The explicit results of these calculations are 
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When accompanied by the auxiliary equation ([21]) for , 
these equations form a closed, autonomous system for the 
orbital elements. 

The results in this section are equivalent to those in 
Sec. Ill C[ which can be easily shown by using Eq. ([39]) 
to replace /* with But they are numerically ill- 

behaved. Specifically, e appears to diverge in the limit 



e ^ 0, and w appears to diverge when sinu = (i.e., 
at every turning point in the orbit). Although these di- 
vergences are canceled analytically by the numerators in 
each case, they are serious obstacles in a numerical inte- 
gration. Thus, the equations given in Sec. Ill Cl are more 
practical, though slightly lengthier. 
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